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On simple cubic lattices, we compute low temperature series expansions for 
the energy, magnetization and susceptibility of the three-state Potts model in 
D = 2 and D = 3 to 45 and 39 excited bonds respectively, and the eight-state 
Potts model in D = 2 to 25 excited bonds. We use a recursive procedure which 
enumerates states explicitly. We analyze the series using Dlog Pade analysis and 
inhomogeneous differential approximants. 
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A subset of the present authors recently described a method [1,2] similar 
to the finite lattice method [3] for generating low temperature series for discrete 
models. This method is based on a recursive computer enumeration of configura- 
tions and has resulted in series expansions for the D = 3 Ising model that extend 
available series by several terms [2-4]. 

In this paper, we present results from a similar analysis for the low temper- 
ature expansions of Potts models in two and three dimensions on a simple cubic 
lattice. We will not describe the method used in much detail. It has already been 
outlined in Ref [2] and will be described in detail in a separate paper. 

The energy for the (/-state Potts model is defined to be 



where <7; is a site-defined field that takes q possible values. The sum is taken over 
all nearest neighbor pairs of spins with 5 being the Kronecker symbol. 

The partition function is the sum of the Boltzmann weights over all config- 
urations 



Sorting configurations by energy, we rewrite this as a sum over E. Defining P(E) 
to be the number of states with a given energy E, we have 

dN 



where d is the number of dimensions, iV is the number of sites and u = e~@ . 

We compute the coefficients P(E) exactly on small systems by recursively 
assembling the system one site at a time. The method enables us to build up 
a lattice with arbitrary length in one direction. Intermediate stages require an 
explicit enumeration of exposed slices transverse to this direction. This effectively 
reduces the computational complexity to that of a system of one less dimension. 

The starting point is a list of all states and corresponding energies for a 
single transverse layer of the lattice. In D = 2, the transverse layer is a line of 
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spins, in D = 3, it is a plane of spins. All the spins outside this transverse layer 
are frozen to the same value; that is, the boundary conditions in the longitudinal 
direction are cold. Spins are then sequentially freed to build up the lattice in the 
longitudinal direction. We store the number of states with a given energy, E, 
and exposed top layer in an array po(E, I), where the integer / is an index which 
specifies the exact configuration of the exposed transverse layer using bit-coding. 
When a new spin is added, we obtain the new counts p' (E, I) as a sum over the 
old counts 

p' (e,i) = J2po(e-a(i,i'),i'). 

I' 

Here I' can differ from / only in the bits representing the newly covered spins, 
and A(J, /') is the change in energy due to any newly changed bonds. For the 
present analysis we add the spins one at a time. Thus, the sum in the above 
equation is only over q terms, representing the q possible values of the newly 
covered spin. After the lattice is grown, a sum over the top layers gives the 
resulting P(E) = J2iPo(E-i -0- We always continue this recursion sufficiently to 
avoid finite size errors in the longitudinal direction. 

As the temperature goes to zero, so does the variable u. Thus, what we have 
in Eq. (3) is the low temperature expansion for Z. From it, we compute the series 
for the average energy, (E) = (u-^) log(Z). Subtracting this expectation value 
before adding the last spin from its value after adding the last spin, we obtain 
the average energy per new site. This also eliminates the effect of the fixed end 
boundaries. Writing, 

(E/N) = ^ 

3 

the low temperature expansion amounts to listing the coefficients ej. 

The recursive technique can be extended to enable calculation of quantities 
such as the magnetization and susceptibility. We define a magnetization in the 
Potts model by 



assigning to each unexcited spin the value one, and to each excited spin the 
value zero. The calculation of susceptibility is carried out using the fluctuation- 
dissipation theorem and we define the low temperature series coefficients Xj as 
follows: 

n x = (m 2 )-(m) 2 = nJ2xjU j 

j 

Let p{E, M, 7) to be the number of states with given energy, magnetization and 
exposed top layer 7. To compute any moment of the magnetization, it would be 
sufficient to compute p(E, M, I). However, one can avoid computing this memory 
expensive quantity. Let us demonstrate this for the case of the magnetization. 

First, note that po(E,I) = ^ M p(E, M, 7) is the count we had be- 
fore. To compute the magnetization, we need one more count: pi(E,I) = 
J2m Mp(E, M, 7). This is because the expectation value of magnetization can be 
written as, 

with P\(E) = ^2jPi (E,I) and Z = J2 E T p (E, I)u E . The counting scheme for 
Pi(E, I) is easy to derive. In analogy with Eq. (4) one can write, 

p[ (E, I) = J2 M P'( E > M ^)=Y1 M ^ E - A e,M- A m , /') = 

M M,I' 

= ( M - A - + A -) P( E - A e , M - A m , 70 = 

M,I' 

= [Pi( E - A e , I') + A m p (E - A e , 70 ] 
r 

Here A e = A e (7, 70 and A m = A m (7, 70 denote the change in energy and 
magnetization when adding the new spin. Thus, computation of the magneti- 
zation series requires just the introduction of one additional count, (which only 
doubles the memory requirement) and we can calculate the magnetization series 
to essentially the same order as the energy series. 

For the susceptibility series, we need to compute (M 2 ). This requires a 
count P2(E, I) = YIm M 2 p(E, M, 7). It is easy to see that p 2 obeys the recursion 



relation, 

p' 2 (E, I) = J2[P2(E-A e , I') + 2A mPl (E - A e , /') + A m 2 Po (E - A e , /') ] (10) 

As discussed in Ref [2], we work on generalized helical lattices and label 
our lattice points by their ordinal number on a helix. In three dimensions, the 
nearest neighbors on the lattice in the x, y and z directions are separated by 
/i x , h y and h z steps along the helix respectively. We assume that the /i's are 
ordered so that h x < h y < h z . Then, our numerical method requires us to keep 
track of at most q hz states and so we try to make h z as small as possible. Let n 
be the effective lattice size, defined as the length of the shortest closed path on 
the helical lattice. For a given set of h values, if we compute the set of nonzero 
vectors S = {n x , n y , n z ; n x h x +n y h y +n z h z = 0} then n = Mms(\n x \ + \n y \ + \n z \) . 
The series expansion will be correct up to the order -u*- 4 " -1 -*. Higher orders are 
corrupted by contributions from graphs that wrap around the lattice. However, 
as described in Ref [2], we can combine results from different helical lattices to 
cancel these finite size effects to some order in the series. In two dimensions, there 
is not enough complexity for this cancellation mechanism to work. Instead, one 
observes that keeping h y spins in the top layer, the optimal choice of the lattice 
is h x = h y — 1. This gives the series correct to order Ah y — 3. 

Our series are listed in Tables I-III. The series for D = 2 and D = 3 
Potts Models were computed on a CM-200/CM-2 Connection Machine using CM- 
Fortran and C* programs. The D = 2, 8-states model series were computed on a 
CRAY-2 using a C code and checked on a CM-2 using CM-Fortran code. To get 
3d series up to 39 excited bonds, we used lattices of effective size up to 10. This 
required the top layer to have at most 15 spins. In Table IV we show the lattices 
and combination factors used. 

Note that our definition of M in Eq. 6 is such that in the completely disor- 
dered state it has the value N/q. The proper order parameter for Potts models 
is the so called reduced magnetization M R which is related to M by the formula 
Mr = [qM — N]/(q — 1). The reduced magnetization takes the values N and 
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in the completely ordered and disordered states respectively. The results we 
give below from our analysis of series are for the reduced magnetization and the 
corresponding susceptibility. 

In addition to the usual dlog Pade (DIP) method [5,6], we will use the method 
of inhomogeneous differential approximants (IDA) introduced by Fisher and Au- 
Yang [7] (see also [8]). These are useful in handling singularities of the form, 

F(u) = A(u)(l -u/u c ) c + B(u) (11) 

where A and B are analytic in u. 

Given a series expansion for F(u) to N-th order, Fn(u) = 1 + J2iLi fi u ' \ 
(we will use the simplification that one can always normalize the series so that 
the constant term is unity), one computes coefficients for polynomials Ql(u) = 
Yh=q Qi u \ Rm(u) = 1 + Yh=i r i u% and Sj(u) = J2i=o s i u \ which satisfy, 

FnQl + Sj = F n Rm (12) 

to order N, with L + M + J = N — 2. Note that for Sj = one gets the usual 
Dlog Pade ratio from Ql/Rm- It is easy to see that potential critical points 
u c are the zeros of Km and for each of these, the exponent £ is estimated as 
C = -Ql{u c )/R' m {u c ). 

Consider first the D = 2 Potts models. Here, we know from self-duality that 
the critical point is at u c = l/(y/q + l). For q < 4 the transition is continuous and 
the critical exponents are known exactly (see [9] and references therein). Models 
with q > 4 undergo a first order phase transition. Having results from both of 
the above categories available, our D = 2 series offer themselves as a good testing 
ground for series analysis methods. 

Given the low temperature series, does one has enough information to de- 
termine the nature of the transition, assuming that the critical temperature is 
exactly known. In D = 2, because of self duality, this is easy if the series at hand 
has a sufficient number of terms. To illustrate this, we plot in Fig. la the energy 
as a function of u from the low temperature series and its dual high temperature 
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series for q = 3 and q = 8. In Fig. lb and Fig. lc, we plot the latent heat L(n) 
derived using duality at the known critical point as a function of the number of 
terms n in the series. The fits of L to a power law in 1/n (Fig. lb and lc) con- 
vincingly demonstrate that the q = 3 model has a second order transition while 
the q = 8 model has a first order transition with the latent heat equal to 1/2 to 
2 parts in a thousand. 

In general however, self duality is not available as a symmetry. In this case, 
one must rely on DIP and IDA analysis on the low temperature series to determine 
the critical properties. Our arguments below are similar in spirit to the discussion 
presented by Enting and Guttmann [10]. 

If the system undergoes a second order phase transition, one expects in gen- 
eral that the order parameter Mr/N vanishes at the critical point, approaching 
it with infinite slope. Estimates of the critical temperature (poles) from D1P- 
s should then cluster well around the exact value and estimates of the critical 
exponent f3 (residues) should also be quite accurate. On the other hand, at a 
first order transition, the magnetization is finite and nonzero and its slope can 
be either finite or infinite. In this case one would expect the approximants to 
continue the curve beyond the critical temperature along the so called pseudo- 
spinodal line [11]. This line intersects the temperature axis at the point us with 
corresponding exponent (3s- Applying DIP-s in this case should then result in a 
systematic overestimation of the critical temperature because it is us that the 
Pade is trying to fit. 

In case of a first order transition with a divergent slope of the magnetization 
as the transition is approached, DIP-s still tend to overestimate the transition 
point because the finite value of the magnetization is not modelled in the DIP-s 
(more detailed reasons can be found in [12]). However, for this case, the IDA-s 
should treat the situation better because they can account for a finite < M/N > 
at the critical point. Thus, comparing the results of the two types of approximants 
one might be able to determine the order of the transition. 

Applying DIP-s to the 45-term magnetization series of the 3-state model in 
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D = 2 leads to a slight systematic underestimation of the critical point. Taking 
into account seven most central approximants we got u c = 0.36595 ± 0.00003 
which is to be compared with the exact value u c = 0.36602.... The error here 
corresponds to the scattering of values from the different DIP-s. In the light of 
the above discussion, this suggests that the transition is continuous. We estimate 
the critical exponent (3 = 0.1084 ± 0.0002 by evaluating it at the known critical 
point for this model. The error bar is of course meaningless as it comes only from 
the error on the extrapolation and ignores the systematic effects of the finiteness 
of the series. The value obtained is about 2.5% below the exact result (3 = 1/9. 

In the 8-state model in D = 2 on the other hand, DIP-s show a critical point 
at u c = 0.2628±0.0003 which is substantially beyond the true value u c = 0.2612.... 
This suggests a first order phase transition. In Fig. 2 we plot u c versus /? for small 
values of J. The points for different J lie fairly well on a line with an obvious 
tendency to overestimate the critical point again. This again establishes the first 
order nature of the transition. The corresponding pseudo-exponent estimated 
from DIP-s has the value (3 S = 0.059 ± 0.005. 

Similar ideas can be applied to the energy and specific heat series. At a 
first order phase transition there is a finite latent heat but the energy curve can 
have either finite or infinite slope (specific heat) as that point is approached. DIP 
analysis of the q = 3 specific heat series in D = 2 shows a slight overestimate 
of the critical point, namely u c = 0.36626 ± 0.00001. IDA-s on the other hand 
lead to a small underestimate (see Fig. 3) giving an overall consistency with 
the second order phase transition present. DIP-s average for critical exponent 
a = 0.412 ±0.001 is rather poor when compared to the exact value a = 1/3. This 
is probably due to the strong confluent singularity present in this case [13]. The 
results of the IDA analysis is shown in Fig. 3 where we plot u c versus a for various 
J values from — 20 with L and M chosen to be equal or differing by at most 
one (see Eq. 12). Notice that if we fit the data to a straight line and compute 
the value of a at the exactly known critical point (vertical line in Fig. 3), we get 
a result which differs from the exact value by about 1%. 
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In the q = 8 model the results from the specific heat series and magnetization 
series are very consistent with each other. There is an overestimate of the critical 
point by DIP-s (us = 0.2620 ±0.0001) as well as by IDA-s. The averaged pseudo- 
exponent from DIP-s is as = 0.592 ± 0.004. 

Finally, an analysis of the susceptibility series for the q = 3 model using 
the Dlog Pade and IDA analysis gave 7 = 1.47 ± 0.02 by extrapolating to the 
known critical point as was done above for a and (3. This is to be compared 
with the exact result 7 = 13/9 = 1.444.... For the q = 8 model, we estimate 
u s = 0.2629 ± 0.0009, 75 = 1.16 ± 0.07. 

Let us now turn to the series for the q = 3 Potts model in D = 3 given 
in Table III. Theoretically, this is the most interesting case of those considered 
in this paper, because of its connection to the SU(3) lattice gauge theory in 
D = 4 [14] and because of the lack of any exact results. There was a good deal 
of confusion about the nature of the transition in the past but by now the first 
order nature of this transition seems to be well established [15]. Although the 
transition temperature is not known exactly, there are very accurate Monte Carlo 
estimates for it. For the purpose of our analysis we will assume that the value 
u c = 0.57659(1) estimated in Ref. [15] is the exact result. We will do so because 
we found that neither the DIP nor the IDA analyses can yield a more accurate 
value. 

Consider first the magnetization series. In Fig. 4 we show the results from 
central Dlog Pades. The data clusters well around the value u c = 0.5785 ±0.0003, 
quite far from 0.57659. IDA-s show the same tendency as can be seen in Fig. 5. 
Here, the results from small J fall very nicely on a straight line beyond the 
critical point which is marked by a cross. These results support the conclusion 
that this model has a first order phase transition in agreement with [10] and 
Monte Carlo data [15]. The critical pseudo-exponent from DIP-s has the value 
Ps = 0.204 ± 0.002 which agrees very accurately with results of Miyashita et al. 
[16] who analyzed a shorter series, and also with numerical simulations [17]. 

Next consider the specific heat series. Here one gets stable results from many 

9 



central Dlog Pades. Also, the IDA-s are quite stable for small J. Fig. 6 shows 
the results of these analyses. The circles correspond to the DIP-s and the other 
symbols are the results from the IDA-s for J < 4. There is no clear evidence 
for systematic overestimation of the critical point by neither DIP-s nor IDA-s 
suggesting that the transition is weakly first order in this variable. The straight 
lines in Fig. 6 are least square fits to IDA-s and DIP-s. Since the latent heat is 
small, one would expect that these should intersect at the critical point where 
they are both dominated by the singularity. Away from the critical point, the 
Dlog Pade and the IDA-s treat the non-leading corrections differently and so 
the results from them could be different. Indeed the lines in Fig. 6 intersect at 
u c = 0.5766(2), a = 0.421(2). We have estimated the error on these parameters 
from the errors in the fitted parameters for the straight lines. 

Finally, we analyzed the q = 3 susceptibility series in three dimensions. 
Here the combined data for DIP and IDA fall nicely on a line. We estimate 
7 = 1.085 ± 0.005 by evaluating the fitted line at u c = 0.57659. 

Recently, Vohwinkel [18] has extended the shadow lattice method and shown 
how one can obtain extremely high order low temperature expansions. His series 
for the magnetization has several more terms than ours and although he does 
not generate series for the other quantities we measure in the present paper, we 
presume he can do so. A challenge now is to see if the ideas of Ref. [18] can be 
incorporated into our method. 
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Figure Captions 

Figure la. The average energy from the series expansions in D = 2 for q = 3 and q = 8. 

Duality was used to get the series in the high temperature phase from the 

series in the low temperature phase. The exactly known transition points 

are shown as vertical lines. 
Figure lb. The Latent heat L(n) as a function of n for q = 3 in D = 2. The solid line 

is a fit to a power law and demonstartes that for n = oo, the latent heat 

vanishes. 

Figure lc. The Latent heat L(n) as a function of n for q = 8 in D = 2. The solid line 

is a fit to a power law plus a constant and demonstartes that for n = oo, the 

latent heat is about 1/2. 
Figure 2. u c versus the exponent (3 from the magnetization series for the q = 8 model 

in D = 2. The exact value of u c is the vertical line. 
Figure 3. u c versus the exponent a from the series for the specific heat for q = 3 in 

D = 2 from IDA analysis. The vertical line is the exact value of u c . 
Figure 4. u c versus the exponent (3 from the magnetization series for the q = 3 model 

in D = 3 using DIP-s. 
Figure 5. u c versus the exponent j3 from the magnetization series for the q = 3 model 

in D = 3 using IDA-s with small J values. The 'exact' value of u c is marked 

with a plus and is a Monte Carlo result from Ref. 15. 
Figure 6. u c versus the exponent a from the series for the specific heat for q = 3 

in D = 3 from DIP and IDA analysis. The transition point is accurately 

determined by the crossing of the lines for DIP-s and IDA-s. 
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Table I: The low temperature expansion coefficients ej, TMj and Xi f° r tne ener gy) magne- 
tization and susceptibility series for the q = 3 Potts model in D = 2 on a simple cubic 
lattice. 



1 














1 





1 











2 











3 











4 


8 


-2 


2 


5 











6 


24 


-8 


16 


7 


28 


-8 


16 


8 


32 


-24 


100 


9 


216 


-72 


216 


10 


160 


-140 


844 


11 


660 


-320 


1,552 


12 


2,072 


-1,164 


7,844 


13 


1,664 


-1,560 


12,112 


14 


11,760 


-7,044 


60,268 


15 


17,700 


-13,000 


118,944 


16 


41,088 


-35,984 


424,072 


17 


156,468 


-101,736 


1,081,392 


18 


207,240 


-219,616 


3,201,728 


19 


849,300 


-647,536 


8,670,688 


20 


1,817,048 


-1,602,194 


25,713,154 


21 


4,021,780 


-3,970,384 


67,206,560 


22 


13,178,264 


-11,239,056 


203,077,760 


23 


25,754,296 


-26,891,584 


532,881,432 


24 


75,653,408 


-73,534,214 


1,558,159,918 


25 


193,458,400 


-191,374,464 


4,250,639,632 


26 


440,725,376 


-486,815,472 


11,956,293,152 


27 


1,296,485,460 


-1,323,802,480 


33,296,697,848 


28 


3,009,317,200 


-3,380,001,144 


92,820,406,096 


29 


7,977,739,920 


-8,964,296,480 


257,249,275,776 


30 


21,217,637,824 


-23,766,809,488 


721,023,458,656 


31 


51,359,965,976 


-61,628,612,552 


1,986,080,278,600 


32 


140,885,970,816 


-165,028,619,666 


5,561,045,323,298 


33 


354,038,121,756 


-432,231,505,864 


15,359,165,767,512 


34 


916,153,258,448 


-1,142,608,252,368 


42,717,426,328,784 


35 


2,439,917,838,708 


-3,039,729,276,192 


118,457,421,095,792 


36 


f -\ f -\ c\r\r\ f \ o a o r\r\ 

6,161,990,034,800 


-7,994,207,679,356 


328,170,466,563,836 


37 


16,397,314,674,708 


-21,295,402,476,752 


909,829,346,983,664 


38 


42,540,620,667,584 


-56,399,959,949,412 


2,520,622,606,225,868 


39 


110,314,458,936,968 


-149,510,058,508,096 


6,973,368,153,491,880 


40 


292,427,669,006,272 


-398,341,255,729,746 


19,322,697,243,220,158 


41 


756,553,239,055,504 


-1,056,154,269,407,136 


53,409,977,638,363,032 


42 


1,994,873,374,110,312 


-2,813,530,068,950,904 




43 


5,238,354,130,103,568 


-7,489,714,245,193,504 




44 


13,686,401,970,717,088 


-19,928,407,714,223,232 




45 


36,195,015,152,016,276 


-53,175,417,534,052,136 
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Table II: The low temperature expansion coefficients €i, TYli and )(i f° r the energy, 
magnetization and susceptibility series for the q — 8 Potts model in D = 2 on a simple 
cubic lattice. 



i 




m, 










1 





1 











2 











3 











4 


28 


-7 


7 


5 











6 


84 


-28 


56 


7 


588 


-168 


336 


8 


-588 


91 





9 


4,536 


-1,512 


4,536 


10 


11,760 


-4,060 


14,504 


11 


-13,860 





15,792 


12 


205,072 


-68,859 


288,169 


13 


144,144 


-84,840 


556,752 


14 


271,460 


-256,424 


2,062,088 


15 


7,553,700 


-2,678,760 


15,132,264 


16 


-713,692 


-2,049,229 


25,582,802 


17 


45,219,048 


-21,023,016 


165,495,792 


18 


232,853,880 


-93,466,856 


720,185,368 


19 


-14,850,780 


-107,162,496 


1,588,846,728 


20 


2,822,644,748 


-1,187,630,969 


10,588,862,669 


21 


6,212,314,080 


-3,159,741,984 


33,856,668,720 


22 


8,166,041,884 


-7,756,117,236 


108,773,186,200 


23 


131,708,763,816 


-56,277,329,304 


596,266,427,232 


24 


167,481,870,528 


-118,516,443,339 


1,709,093,729,238 


25 


846,878,642,400 


-506,752,816,584 


7,126,592,218,032 



15 



Table III: The low temperature expansion coefficients ej, Tfli and Xi f° r the energy, 
magnetization and susceptibility series for the q — 3 Potts model in D = 3 on a simple 
cubic lattice. 



i 




m, 










1 





1 











2 











3 











4 











5 











6 


12 


-2 


2 


7 











8 











9 











fO 


60 


-12 


24 


ff 


66 


-12 


24 


12 


-168 


28 


-56 


13 











14 


420 


-90 


270 


15 


900 


-180 


540 


16 


-1,728 


318 


-930 


17 


-2,448 


432 


-1,296 


18 


6,708 


-1,320 


4,768 


19 


9,462 


-1,992 


7,968 


20 


-14,280 


2,760 


-10,560 


21 


-49,686 


9,368 


-36,992 


22 


71,940 


-14,460 


64,812 


23 


177,192 


-35,280 


163,440. 


24 


-194,544 


36,680 


-16,5464 


25 


-684,300 


134,568 


-659,088 


26 


515,892 


-108,516 


600,024 


27 


3,087,234 


-609,692 


3,278,256 


28 


-1,927,296 


370,500 


-1,980,408 


29 


-10,943,904 


2,153,016 


-12,285,816 


30 


3,863,712 


-792,218 


5,005,014 


31 


44,383,506 


-8,867,580 


55,200,864 


32 


-4,406,976 


935,124 


-6,062,712 


33 


-177,069,948 


34,889,512 


-227,203,096 


34 


-1,133,220 


63,834 


1,954,650 


35 


652,560,090 


-130,265,472 


914,339,736 


36 


199,263,288 


-39,322,372 




37 


-2,553,456,210 


507,892,056 




38 


-1,235,636,652 


239,776,590 




39 


9,742,992,324 


-1,940,344,524 
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Table IV: The lattices parameters and combination factors that give the series accurate 
to 39 excited bonds in D = 3. 



h x 


hy 


h z 


Coefficient 


9 


14 


15 


2 


11 


12 


15 


-1 


9 


11 


15 


-2 


10 


13 


14 


1 


11 


12 


14 


5 


9 


11 


14 


-1 


7 


12 


13 


1 


10 


11 


13 


-3 


8 


10 


13 


1 


5 


11 


12 


3 


7 


10 


12 


-5 
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